%% 1. Estimate parameter values for different higher order risk attitudes

clear

% Load data
T = readtable('data_reflection_horp.xlsx');

% Uncomment line below to perform analysis for original data only
T(123:245,:) = [];

% Calculate number of risk apportionment choices per subject
apport_choices = zeros(size(T,1),9);
treatments = {'loss','gain','sprg'};
for i=1:3
for j=1:3
for k='a':'l'
ind = strcat(treatments{i},int2str(j+1),'_',k);
apport_choices(:,i*3-3+j) = apport_choices(:,i*3-3+j) + T.(ind);
end
end
end
% The first three columns of apport_choices contain the number of risk
% apportionment choices for order 2, 3 and 4 for losses, followed by the
% same for 50-50 gains and small probability gains.

% Aggregate responses
F = [];
for i=1:9
F = [F; histcounts(apport_choices(:,i),13)];
end
F = flip(F,2);

no_msg = optimset('Display','off'); % disable messages about convergence

% Set up initial values and constraints
x0 = [0,0,0,0,0];
Aeq = [1,1,1,0,0,;0,0,0,0,0;0,0,0,0,0;0,0,0,0,0;0,0,0,0,0];
Beq = [1,0,0,0,0];
lb = [0,0,0,0,0];
ub = [1,1,1,0.5,0.5];
r_Aeq = [1,1,1,0,0,;0,0,0,0,0;0,0,0,0,0;0,0,0,0,0;1,0,-1,0,0];

for i = 1:3
    for j = 0:2
	% Estimate maximum likelihood values and unrestricted likelihood
	[B(i*3-2+j,:),l] = fmincon(@(v)mle_overall(v,F(i+j*3,:)),x0,[],[],Aeq,Beq,lb,ub,[],no_msg);
	u = -l;

	% Calculate the restricted likelihood 
	[c,l] = fmincon(@(v)mle_overall(v,F(i+j*3,:)),x0,[],[],r_Aeq,Beq,lb,ub,[],no_msg);
	r = -l;

	% Calculate the log likelihood ratio.
	LR(i*3-2+j,1) = -2*(r - u);
    LR(i*3-2+j,2) = 1-chi2cdf(-2*(r - u),1);
    end
end

disp('Maximum Likelihood Estimations:')
disp(B)
disp('Likelihood Ratios and p-values, with restriction lambda_s = lambda_o:')
disp(LR)

%% 2. Test reflection effect

clearvars -except apport_choices F

no_msg = optimset('Display','off');

x0 = [0,0,0,0,0,0,0,0,0,0];
Aeq = [1,1,1,0,0,0,0,0,0,0;0,0,0,0,0,1,1,1,0,0;zeros(8,10)];
Beq = [1,1,0,0,0,0,0,0,0,0];
lb = [0,0,0,0,0,0,0,0,0,0];
ub = [1,1,1,0.5,0.5,1,1,1,0.5,0.5];

k=1;
for i = 1:3
    for j = [3, 6]
	% Estimate maximum likelihood values and unrestricted likelihood
	[B(k,:),l] = fmincon(@(v)mle_reflection(v,F(i,:),F(i+j,:)),x0,[],[],Aeq,Beq,lb,ub,[],no_msg);
	u = -l;

	% Calculate the restricted likelihood 
	[C(k,:),l] = fmincon(@(v)mle_reflection_rsc(v,F(i,:),F(i+j,:)),x0,[],[],Aeq,Beq,lb,ub,@mle_reflection_con,no_msg);
	r = -l*10^-5; %mle_reflection_rsc is rescaled for convergence

	% Calculate the log likelihood ratio.
	LR(k,1) = -2*(r - u);
    LR(k,2) = 1-chi2cdf(-2*(r - u),1);
    k=k+1;
    end
end

disp('Likelihood Ratios and p-values, with restriction lambda_s(g)/lambda_o(g) = lambda_s(l)/lambda_o(l):')
disp(LR)

%% 3. Estimate combinations of second and fourth order risk preferences

clearvars -except apport_choices F

% Create arrays without missing observations
for i=1:3
second_order(:,i) = {apport_choices(~isnan(apport_choices(:,i*3-2)),i*3-2)};
end
for i=1:3
fourth_order(:,i) = {apport_choices(~isnan(apport_choices(:,i*3)),i*3)};
end

% Set up initial values and bounds
x0 = [0.125,0.125,0.125,0.125,0.125,0.125,0.125,0.125,0.07,0.07,0.07,0.07];
lb = [0,0,0,0,0,0,0,0,0,0,0,0];
ub = [1,1,1,1,1,1,1,1,0.5,0.5,0.5,0.5];
r_x0 = [0.125,0.125,0.125,0.125,0.125,0.125,0.125,0.07,0.07,0.07,0.07];
r_lb = [0,0,0,0,0,0,0,0,0,0,0];
r_ub = [1,1,1,1,1,1,1,0.5,0.5,0.5,0.5];

for i=1:3
    % Create variable with data to pass to mle function
	d{i} = second_order{i} + fourth_order{i}*100;

	% Estimate maximum likelihood values
	B(i,:) = mle(d{i},'pdf',@mle_mixedrisk,'start',x0,'lowerbound',lb,'upperbound',ub,'optimfun','fmincon','options',statset('FunValCheck','off'));

	% Calculate the unrestricted likelihood.
	estimator_vals{i} = num2cell(B(i,:));
	u(i) = sum(log(mle_mixedrisk(d{i},estimator_vals{i}{:})));

	% Calculate the restricted likelihood.
	mle(d{i},'pdf',@mle_mixedrisk_rsc,'start',r_x0,'lowerbound',r_lb,'upperbound',r_ub,'optimfun','fmincon','options',statset('FunValCheck','off'));
	r_vals{i} = num2cell(ans(:));
	r(i) = sum(log(mle_mixedrisk_rsc(d{i},r_vals{i}{:})));

	LR(i,1) = -2*(r(i) - u(i));
    LR(i,2) = 1-chi2cdf(LR(i,1),1);
end

% Insert the paramater values backed out from the constraint
B = [B(:,1:8),1 - sum(B(:,1:8),2), B(:,9:12)];

disp('Maximum Likelihood Estimations:')
disp(B)
disp('Likelihood Ratios and p-values, with restriction lambda_ss + lambda_oo = lambda_so + lambda_os:')
disp(LR)